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THE  NONUNEAR  DYNAMICS  OF  MHD  UNSTABLE  PLASMAS 
X.  Quasi-Llnear  Theory  of  MHD  Instabilities 

This  chapter  begins  the  second  part  of  this  book,  that  concerning  the 
nonlinear  theory  of  MHD  instabilities.  Before  commencing  it  is  necessary 
to  point  out  that  the  nonlinear  theory  of  MHD  instabilities  is  not  nearly 
as  well  developed  as  the  linear  theory  and  much  of  what  will  be  discussed 
here  is  necessarially  more  speculative  than  the  material  discussed  in  the 
preceeding  chapters,  which  treated  the  linear  theory. 

We  begin  by  examing  the  quasi-linear  theory  of  MHD  instabilities  in 
which  the  perturbed  velocities  are  non-singular.  This  applies  potentially 
to  id  s  1  kink  tearing  modes  and  double  tearing  modes  (Chapters  VIII  E  and 
F)  as  well  as  to  ideal  MHD  modes.  In  this  and  the  next  chapter  we  change 
the  notation  slightly  from  that  in  the  previous  chapters.  Now  any  quantity 
(say  A)  will  be  denoted  by(A)  +  A  where(A)  is  the  ensemble  average  and  A  is 
the  perturbation  about  the  ensemble  average.  Ensemble  averaging  will  be 
denoted  by  triangular  brackets 

For  instance,  if  we  restrict  ourself  to  cylindrical  geometry,  A  is  a 


function  of  radius  only,  and 

rj 


(X  1) 


the  founder  the  summation  indicating  m  and  k  are  not  both  *  0.  Since  any 
physical  quantity  must  be  real,  the  complex  conjugate  is  added  on  in  Eq. 

(X  1).  The  process  of  taking  the  ensemble  average  can  then  be  defined  as 
an  average  over  8  from  zero  to  2ir  and  an  average  over  z  from  0  to  L,  where 
L  is  either  a  periodicity  length  if  such  a  length  exists  or  else  is  some 
very  large  length.  Clearly  (A)  •  0,  but 

<a  5>= rxcr,»,A)  +  c  c 


£ — / 

A, 1**0 


(X  2) 
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Expressing  the  fluid  quantities  in  this  way,  the  four  fluid  equations  for 


i«s>+?)=  i*  (<y>+ m 

&%)'***)+  f(<r>*VS' C<V>  +  ^  *°w 


where  we  have  assumed  y  «  5/3  to  avoid  confusion  between  the  adiabadic  y 
and  the  growth  rate. 

The  first  step  is  to  take  the  ensemble  average  of  (Eqs.  X  3).  The 
result  is 


+  v- py  =  -  v  <j> 

*  (a) 

iff^t  «p>- 

<3y>-v<y>  -  <(5<y>-2?>  ~f> "  <^>  ysv> 


(X  4) 
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(b) 


*£>  -  7x<v>x(?)s  ?* 


+  <v><KP>*  fW-W* 
*  '  " 


(c) 


T 


<?*•*> 


(d) 
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^  •  ■ 

Then  subtracting  Eqs.  (X  4)  from  Eqs.  (X  3)  gives  the  equations 

;  • 

i 

j.  S<v>.v<v>*,oV-7<v>»<»><y>-2v»C? 

■  !  'rtf  '  i/  1  '  “ 

?)*<8> ■*' (?*<&>) *B)y  “? 

.  }  $.  J  <y>  +  <P  y  >.?<v>  -  ^>-5  5  *■  <f  «>•*  -v> 

_y,>y.  *  v  «.  <^><  y  •  h  y>  -  ^y-sy  +  <f'‘--^ 

+  u_(5,|>5- £«?*«*  ?>  ^ 

(X  5)  | 

-  7xVx<B>-  yx<v>x?=  ^  x  <  v  > 

to  j 

|  H  +  V-  7<p>  +  <V>  -7?  +  f  ?  ?-<v>  ♦  £  <  P>  V-  v  » 

|  .  +  <S-V  JO'f  +  #<?  7<5> 
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(*0  j 

Equations  (X  4  a-d)  and  Equations  (X  5  a-d)  are  now  a  complete  set  of  non¬ 
linear  equations  for(p),p, V,  (b),  B,  (p)and  p,  and  so  far  no  approximations 
have  been  made. 

The  quasi-linear  approximation  to  these  equations  consists  of  neglecting 
the  (  pV*VV )  term  (that  is  the  term  cubic  in  fluctuating  quantities)  from 
the  right  hand  side  of  Eq.  (X  4b)  and  in  neglecting  the  entire  right  hand 
side  of  Eqs.  (X  4  a-d).  We  will  briefly  discuss  these  approximations  now. 

The  fluctuating  quantities  are  assumed  to  be  small  so  that  a  term  like 
(pV*7V)  is  expected  to  be  much  smaller  than  a  term  like  p^V'Vv)  since 
p«p.  This  is  the  basic  justification  for  neglecting  the  cubic  nonlinear 
terms  on  the  right  hand  side  of  Eq.  (X  4b). 

Notice  that  all  terms  on  the  right  hand  side  of  Eqs.  (X  5  a-d)  are 
various  quantities  minus  their  ensemble  averages.  This  ensures  that  for 
our  cylindrical  case  of  azimuthal  and  axial  averaging,  there  is  no  quantity 
on  the  right  hand  side  of  Eqs.  (X  5a-d)  which  is  a  function  of  only  r.  That 
is  every  driving  term  on  the  right  hand  side  has  oscilliatory  structure  in 
0  and/or  z.  For  instance  if  there  are  two  fluctuations  at  (k^.m^)  and  ^.n^), 
the  quadratic  terms  on  the  right  hand  side  will  drive  additional  fluctuations 
at(2kp  2m^),  (2k2>  2m^)  and  (k^  ±  k2,  m^  ±  m^) .  However  no  fluctuation  will 
be  driven  at  k  *  a  ■  0  since  the  ensemble  (ie  azimuthal  and  axial)  average 
average  of  the  right  hand  side  is  subtracted  out.  These  terms  then  describe 
the  coupling  of  fluctuations  at  different  wave  lengths  directly  with  each 
other. 

Thus  the  quasi-linear  approximation  neglects  coupling  of  the  fluctuations 
among  themselves,  but  Includes  the  coupling  (correct  to  second  order)  of  the 
fluctuating  to  ensemble  average  quantities.  Therefore  the  basic  idea  behind 
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the  quasi-linear  approximation  is  that  coupling  of  the  waves  to  the  background 
is  somehow  more  important  than  the  coupling  of  the  waves  to  each  other. 


With  the  right  hand  side  of  Eqs.  (X  5  a-d)  set  equal  to  zero,  the 

equations  for  the  tilda'd  quantities  are  very  similar  to  the  linearized  fluid 

equations  of  earlier  chapters,  but  with  two  important  complications.  First 

of  all  the  background  velocity  is  not  equal  to  zero  in  Eqs.  (X  5  a-d)  and 

secondly,  the  background  quantities  are  functions  of  time  as  well  as  radius. 

To  proceed,  we  now  restrict  ourselves  to  the  lowest  order  tokamak  ordering 

and  to  systems  not  far  above  stability  threshold.  As  we  will  now  show,  these 

complications  can  then  be  eliminated .In  lowest  order  tokamak  ordering,  recall 

Bz»^  ,  R»a  and  the  mode  structure  is  two  dimensional  in  the  r8 plane  (that 

is  V  “  B  *  0).  Also  V»V  =0,  that  is,  the  perturbation  is  incompressible, 
z  z 

The  equation  for  B  ,  from  Eq.  (X  4c)  then  becomes 
z 


2£*_  _L  -2-  rv,  8, 

M  -  '  r  9r 


(X  6) 


Since  B  is  very  large,  it  cannot  change  very  much  without  drastically 
z 

altering  the  systems  energy  density.  However  for  B  to  remain  unchanged, 

z 


it  must  be  that 


Vr  Z  O 


(X  7) 

Thus  to  lowest  order  in  tokamak  ordering,  the  ensemble 
average  radial  velocity  vanishes.  Now  let  us  consider  the  azimuthal  and 
axial  component  of  velocity.  Making  use  of  the  fact  that  the  fluctuating 
flow  and  magnetic  field  are  incompressible  in  two  dimensions,  it  is  a 
straight  forward  matter  to  show  that  the  6  and  z  components  of  the  terms 


3V_ 


'8 


on  the  right  hand  side  of  Eq.  (X  4b)  are  proportional  to  ^  and  . 
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Thus  there  is  no  term  which  acts  as  a  source  term  for  an  acceleration  in 
6  and  z  (one  can  also  show  that  this  is  true  in  cylindrical  geometry  without 
an  expansion  in  tokamak  ordering).  Therefore,  if  V.  and  V  are  initially 
zero,  they  will  be  zero  at  all  subsequent  time.  Hence,  to  lowest  order  in 
tokamak  ordering,  there  is  no  induced  fluid  velocity,  or 


V  =  o 


(x  8) 


(In  general  in  cylindrical  geometry  “0,  but  V^_  0.) 

In  this  case  the  equation  for  the  perturbed  quantities  (Eqs.  (X  5  a-d) 
with  the  right  hand  side  m  0)  are  different  from  the  linearized  equations 
in  previous  chapters  only  in  that  the  background  quantities  are  functions 
of  time.  Therefore  the  perturbed  quantities  cannot  be  assumed  to  have  a 
time  dependence  like  exp  yt.  However  if  the  time  scale  for  the  change  of 
equilibrium  quantities,  t  satisfies  the  condition 


i 


(X  9) 


where  y  is  the  linear  growth  rate  (assuming  the  ensemble  average  quantities 
are  time  invariant ) .then  one  can  solve  Eqs.  (X  5  a-d)  by  making  a  WKB 
approximation.  That  is,  fluctuating  quantities  are  proportional  to  exp J y(t)dt. 
The  zero  order  WKB  approximation  to  Eqs.  (X  5  a-d)  then  give  the 
conventional  result  ,  . 

O-N 

if  +  v*  =  o  , 

,  *  (?*««>*»>$  CM 

(‘) 


M>«)> 

A) 

y  B  s 


V 


Vx  vx<BW> 


(X  10) 


7*  V  -  o 


M) 
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where  we  have  assume  the  flow  is  Incompressible  in  Eq.  (X  lOd).  These  are 
just  the  conventional  equations  for  perturbed  quantities  which  we  have 
studied  in  previous  sections.  Using  them  to  express  p,  B  and  p  in  terms 


where  $r  is  the  radial  displacement,  f  ■  Vr/y  and  where  as  usual  the 

perturbation  is  assumed  to  be  a  summation  of  individual  fluctuations  proportional 

to  exp  (im  0  +  ikz) .  It  only  remains  to  check  a  postiori  that  Eq.  (X  9) 

is  satisfied.  If  r  is  the  radial  scale  length  for  say  the  density  variation, 
n 

Eq.  (X  9)  reduces  to 

Cr  <  <  L  n  (X  12) 


Thus  Eq.  (X  9)  is  satisfied  as  long  as  the  radial  displacement  is  small 
compared  to  the  radial  scale  length. 

Equations  (X  //  a-c)  then  describe  the  evolution  of  the  background 
density,  pressure  and  poloidal  field  in  response  to  the  instability.  These 
equations  in  themselves  do  not  guarantee  pressure  balance.  However  pressure 
balance  can  easily  be  maintained  through  very  small  changes  in  Thus 

we  expect  that  will  not  exactly  vanish,  but  rather  that  there  will  be  a 
very  small  which  will  modify  very  slightly,  (according  to  Eq.  (X  6  )) 

so  as  to  maintain  pressure  balance  at  all  times. 


7 


To  summarize,  the  quasi-linear  approximation  is  that  the  instability 
grows  at  its  temporarilly  local  linear  growth  rate  and  that  the  background 


plasma  slowly  reacts  to  the  presence  of  this  instability,  modifying  its 
growth  rate.  To  our  knowledge,  no  general  conditions  under  which  quasi- 
linear  theory  is  valid  have  been  given.  Clearly  it  is  not  valid  in  a 
situation  like  fluid  turbulence  where  the  main  interaction  is  the  modes 
with  each  other  rather  than  the  modes  with  the  background.  The  quasi-linear 
theory  is  most  likely  to  be  valid  in  a  situation  where  a  plasma  is  initially 
weakly  unstable  and  where  the  subsequent  growth  of  instability  causes  it 
to  evolve,  according  to  Eqs.  (X  11  a-c)  to  a  stable  state.  It  is  least 
likely  to  be  valid  in  a  condition  of  steady  state  or  fully  developed 
turbulence. 

An  interesting  related  observation  is  that  the  quasi-linear  theory  of 
velocity  space  instabilities  has  been  widely  studied  for  infinite  homogeneous 
plasmas.  There,  comparisons  between  quasi-linear  theory  and  particle 
simulation  has  frequently  shown  that  quasi-linear  theory  does  provide  a 
reasonably  accurate  approximation  of  the  description  of  the  evolution  of 
an  initially  unstable  plasma  to  a  final  stable  state. 

There  is  one  aspect  of  Eqs.  (X  11  a  and  c)  which  might  at  first  appear 
paradoxical.  That  is  the  ensemble  average  magnetic  field  is  no  longer 
frozen  into  the  ensemble  average  flow.  This  can  most  easily  be  seen  in 
slab  geometry.  If  the  perturbation  is  two  dimensional  and  incompressible 
in  the  xy  plane,  a  simple  calculation  shows  that  the  slab  geometry  analogs  to 
Eqs.  (X  11  a  .  and  c)  are 


a  r<£> 

»  a* 

y 

PA4 
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o) 

(0 


(X  13) 


Since  Vx  =  0,  frozen  in  field  means  p/By  is  constant  which  clearly  is  not 
the  case  if  p  and  By  obey  Eqs.  (X  13). 

The  resolution  of  this  apparent  paradox  comes  from  the  fact  that  the 

MHD  constraint  does  not  require  the  average  field  to  be  frozen  into  the 

average  density;  rather  it  requires  that  the  exact  field  be  frozen  into 

the  exact  flow.  The  two  requirements  are  not  the  same  as  we  will  now 

show.  To  show  this  imagine  a  fluid  with  initial  density  PQ(x)  and  field 

B  (x)  iv*  Then  give  each  element  of  the  fluid  a  displacement  £  sin  ky  i 
o  -y  ~x 

where  £  is  constant.  Let  us  now  calculate  the  new  density  and  magnetic 
field.  To  do  consider  3  neighboring  points  initially  at  (x,y),  (x  +  dx,  y) 
and  (x,  y  +  dy).  After  displacement,  these  points  go  to  (x  +  £  sin  ky,  y), 
(x  +  dx  +  £  sin  ky,y),  and  (x  +  £  sin  ky  +  kdy  cos  ky,  y  +  dy) .  The  two 
difference  vectors  before  and  after  the  displacement  are  dxix  and  dy  iy 
before;  and  dxl  and  £  dy  k  cos  ky  i  +  dy  iy  after.  The  original  and 
final  area  is  dxdy,  as  can  be  seen  by  taking  cross  products  of  the  two  sets 
of  difference  vectors.  This  confirms  the  incompressible  character  of  the 
displacement.  Thus  if  the  initial  position  of  a  point  (ie  (x»y))  is 
denoted  rQ,  and  the  final  position  (x  + £  sin  ky,y)  is  denoted  r,  then 


(X  14) 


Solving  for  rQ  in  terms  of  r,  we  have 


p(r)  =  ^  (x  "  **"^3  ) 


(X  15) 


Assuming  that  the  variation  of  pQ  is  small  in  a  distance  {,  we  find  that 
the  density  averaged  over  y  is  given  approximately  by 

s  .  \  -  r..\  .  I  ff^ 


<pw>  =  i  i  $ 


(X  16) 
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We  will  now  continue  by  calculating  the  field,  exploiting  the  fact 
that  it  is  frozen  into  the  flow.  In  the  undisplaced  fluid,  the  field  is 


parallel  to  the  first  difference  vector  dy  i^.  This  difference  vector  is 
displaced  to  dy  (i  k  cos  ky  i  +  iy)  and  this  must  be  parallel  to  the 
'displaced'  field.  Thus  a  unit  vector  parallel  to  B  is 

f  Jk  raa>  i  x  +  i  3 
'8~  (h  t'j?  tm?  )"■ 


(X  17) 


Now  calculate  the  magnitude  of  B.  The  initial  flux  through  the  line 


{(x+  dx,  y)  -  (x,y)>  is  (B  (x)i  )*(dxi  )  -  B  (x  )dx.  After  displacement 

o  — y  — y  o  o 

flux  is  frozen  in,  so 

8()<8-  l)M  "  or  , 


B(r)  Be1-*' 


)((fj 


The  y  component  of  B  in  the  displaced  fluid  is 

B  =  8(0  iV  ^ 

Again,  assuming  that  Bq  varies  only  slightlyin  a  distance  £,  we  find 


(X  19) 


<8„>=  8.00  + 


(X  20) 


As  is  apparent  from  Eqs.  (X  16  and  20),  ^p)/^By)  t  n/B  so  the  average  field 
is  not  frozen  into  the  average  density.  Yet  the  field  was  calculated  by 
making  explicit  use  of  the  fact  that  the  exact  field  is  frozen  into  the 


exact  density.  Therefore  ensemble  averaging  destroys  the  'frozen  in'  nature 
of  magnetic  fields  in  ideal  MHD. 
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In  a  sense,  this  is  not  a  surprising  result.  As  we  have  seen  in 


i 


Chapter  VII  and  VIII,  fluid  flow  with  frozen  in  field  can  lead  to  very 
complicated,  small  scale  magnetic  structures  which  may  well  be  completely 
washed  out  with  any  kind  of  ensemble  averaging.  However,  this  complicated 
structure  was  a  direct  consequence  of  frozen  in  fields,  so  averaging  field 
and  density  will  in  all  likelyhood  destroy  this  link. 

An  analogous  situation  is  an  incompressible  fluid  with  variable  density 
p(x)  filling  the  square  -L  <  x.y.  <  L.  Imagine  stirring  up  this  fluid  for  a 
long  time.  Even  though  the  flow  is  imcompressible,low  and  high  density 
parts  are  forced  near  each  other  in  an  almost  random  way.  Hence  any  kind 
of  coarse  grain  averaging  of  the  fluid  will  give  rise  to  nearly  uniform  density, 
even  though  this  might  at  first  sight  appear  to  be  impossible  due  to  the  incom¬ 
pressibility  condition.  Mixing  scotch  and  soda  is  a  more  familiar  example  of 

this  effect.  Therefore,  in  all  cases,  ’microscopic'  constraints  can  become 
unglued  upon  ensemble  averaging. 

Let  us  now  apply  the  quasi-linear  theory  to  an  internal  m  *  1  kink 
tearing  mode.  The  theory  is  greatly  simplified  because  *  constant  (that 
is  independent  of  r)  within  the  q  -  1  surface,  and  only  the  m  =  n  *  1  mode 
is  assumed  to  exist.  As  it  can  easily  be  shown  that  for  B  »  constant  and 
Cr  “  constant,  there  is  no  contribution  to  from  the  equation  for 


<06 


e 


(X  21) 


which  is  very  much  like  a  diffusion  equation.  Another  way  of  writing 
Eq.  (X  21)  is  to  take  the  curl  and  write  it  as  an  equation  for  the  plasma 
current  in  the  Z  direction.  The  result 


(X  22) 


i 

] 
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which  Is  a  diffusion  equation  for  Jz. 

We  conclude  by  examing  what  sort  of  an  interpretation  for  m  ■  n  ■  1 

sawtooth  oscillations  can  be  made  by  using  quasi-linear  theory.  Imagine 

that  at  time  t  -  0,  an  Ohmically  heated  tokamak  discharge  has  q(r  ■  0)<  1; 

and,  as  is  usually  the  case,  monotonically  decreasing  profiles  of  p(r),  p(r) 

and  Jz(r)  (but  monotinically  increasing  profile  of  q(r))  as  shown  in  Fig.  (X  la). 

This  configuration  is  unstable  to  the  m  *  n  =  1  kink  tearing  mode.  According 

to  quasi-linear  theory,  this  mode  will  grow  at  its  local  growth  rate  and  p, 

p  and  J  all  obey  diffusion  equations.  Thus  current,  density  and  pressure 
z 

all  diffuse  outward,  and  tend  to  evolve  the  plasma  to  a  final  state  where 
Jz  (and  thereby  q),  p  and  p  all  have  no  radial  gradient  within  the  singular 
region,  as  shown  in  Fig.  (X  lb).  Now  let  us  postulate  that  the  steep  gradient 
near  the  singular  surface  in  Fig.  (X  lb)  rapidly  diffuses  away,  either  by 
classical  diffusion,  or  more  likely  by  exciting  some  sort  of  micro-instability. 
After  this,  the  profile  has  smoothed  out  and  looks  like  that  shown  in 
Fig.  (X  lc),  where  the  interaction  between  the  m  =  n  =*  1  instability  and 
whatever  smoothes  out  the  profile,  will  insure  that  in  the  final  state, 
q(r  >■  0)^1,  so  the  plasma  is  stable.  However  the  plasma  is  heated  by  Ohmic 
heating,  so  that  the  hotter  central  region  tends  to  be  preferentially  heated. 

This  lowers  the  resistivity,  and  channels  the  current  into  the  center 
thereby  lowering  q(r  ■  0)  until  it  is  less  than  unity.  Thus  the  current 
channeling  tends  to  force  the  plasma  back  into  an  unstable  state.  This 
channeling  plus  the  quasi-linear  reaction  of  the  plasma  are  the  two  basic 
components  of  a  relaxation  oscillation,  very  much  as  observed  in  the  central 
region  of  tokamaks. 
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Let  us  now  see  how  large  would  be  for  such  a  relaxation  oscillation. 
From  (Eq.  (X  lib) the  energy  confinement  time  of  the  central  region  is  given 

-  -I 

roughly  by  p  i 

«-  •'(f) 

'  V  J  (X  23) 

where  rs  is  the  radius  of  the  singular  surface.  Let  us  consider  for  instance 

TFR  where  Tg  -  2Kev,  *v»  -j,  rg  s  10  cm  and  F1  ^  2.  Then  according  to 

4  S 

Eq.  (VIII  44),  y  'v*  10  .  The  fast  downstroke  of  the  relaxation  oscillation 

occurs  in  about  one  millisecond.  Even  assuming  all  of  the  energy  is  expelled 

from  the  central  region  in  this  time,  Eq.  (X  lib)  gives  (t  /r  )  'V/  i.  Since 

r  s  j 

much  less  than  all  of  the  energy  is  expelled  (recall  that  the  diffusion  only 

flattens  the  pressure  profile  within  the  singular  surface),  r/r  is  really 

s 

much  less  than  one  third,  so  at  least  this  aspect  of  quasi-linear  theory 
appears  to  be  valid.  In  a  later  section  we  will  discuss  numerical  simulations 
of  this  instability.  While  many  of  the  gross  features  can  be  described  by 
quasi-linear  theory,  we  will  see  that  the  actual  description  can  be  quite  a 


bit  more  complicated. 
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X  1)  The  quaal-llnear  evolution  of  the  toroidal  current  for  in  n  *  n  ■  1 
internal  kink  tearing  mode,  a)  Ac  initial  current,  b)  The  current  on 
completion  of  the  MHD  inatability  phaae,  c)  The  current  efter  the  reault- 
lng  ateep  gradlenta  dlffuae  away. 


si-Llnear  Theory  and  Rutherford ' s  Nonlinear  Theory  of  Tearing  Modes 


Since  the  tearing  node  is  driven  by  magnetic  energy  in  the  outer  region, 

it  is  tempting  to  describe  its  nonlinear  evolution  in  terms  of  a  quasi- linear 

theory  of  only  the  outer  region.  In  this  way,  one  might  be  able  to  follow 

the  growth  of  the  tearing  mode  coupled  to  the  depletion  of  free  energy  which 

drives  it.  Then,  when  all  of  the  free  energy  is  used  up,  linear  growth  should 

stop.  Another  attractive  feature  of  an  outer  region  quasi-linear  theory  is 

that  it  works  the  same  way  no  matter  what  the  dissipation  mechanism  in  the 

inner  region  is.  The  problem  is  that  the  MHD  equation  for  the  outer  region 

is  singular,  so  that  one  still  has  to  treat  the  singularity.  For  instance  if 

B  (x,y,t)  ■  B  (x)  exp  i  (ky  +  yt)  +  c.c.,  then  Eq.  (X  13b)  reduces  to 
*  X  a. 


?<y  Jr  >*  1^1  .  JL  ^  .  £' '  Sf  *  r.  c. 
>/  *  *  ** 


(XI  1 


where  we  have  related  B^  to  the  displacement  by  the  ideal  MHD  relation  B^  ■ 
i  k  (By)  (.  Also,  as  in  Chapter  X,  we  have  assumed  the  presence  of  a  large, 
uniform  B^.  Clearly,  Eq.  (XI  1)  is  singular  at  x  ■  0,  the  position  where 
(By)  ■  0.  For  (By (x))  as  shown  in  Fig.  (XI  la),  it  is  possible  to  show 
that  in  the  outer  region,  the  ambient  magnetic  field  loses  energy  and  thereby 
drives  the  mode.  To  do  so,  make  use  of  the  fact  that  around  x  ■  0,  because 
of  the  nonzero  resistivity  (  (x)  is  actually  a  well  behaved  function  of  x 
over  the  entire  domain  -  ®  <  x  <  ®.  Therefore,  the  last  form  of  Eq.  XI  1 
can  be  Integrated  by  parts  so  that 

$.&<»&*'■  - 

.  n.  a  •/: 7  ;S<V  " 


r 


Since  (By)  ■  0  at  x  ■  0,  the  last  Integral  in  Eq.  X  2  is  not  singular. 

Also  By  and  its  second  derivative  everywhere  have  opposite  signs  (see  fig. 

(X  la  and  b),  so  it  is  clear  that  energy  is  released  from  initial  magnetic 
field.  This  energy  goes  to  drive  the  tearing  mode.  Clearly  the  energy 
released  within  the  inner  region  is  not  accurately  calculated;  however, 
this  should  be  unimportant,  very  little  magnetic  energy  is  there  in  the 
first  place. 

We  now  examine  how  the  quasi-linear  evolution  of  the  plasma  in  the 
outer  region  can  be  carried  out.  There  are  two  relevant  widths  in  the  inner 
region,  first  there  is  the  resistive  layer  size  given  by  Eq.  (VII  25c) 
and  secondly  there  is  the  island  width  Ax  .  It  is  clear  that  Eq.  (XI  1) 

IS 

is  not  valid  for  x  <  L  because  all  information  concerning  resistivity  was 

c 

left  out.  However,  it  is  also  true  that  Eq.  (XI  1)  cannot  be  valid  for 
x  <  either.  The  reason  is  that  it  is  a  partial  differential  equation 

relating  the  change  in  (B^  to  the  spatial  structure  of  (By  )  in  the 
immediate  neighborhood.  However,  during  reconnection,  the  magnetic  field 
line  is  suddenly  affected  by  the  other  field  line  it  reconnects  with;  a 
field  line  which  is  generally  very  far  away.  Thus  the  behavior  of  the 
ambient  field  during  reconnection  cannot  be  described  by  a  single  partial 
differential  equation.  This  can  also  be  seen  in  the  derivation  leading  up 
to  Eq.  (X  20).  From  this  it  is  clear  that  Eq.  (XI  1)  described  the  response 
of  the  average  frozen  in  field  to  small  displacements  of  the  fluid. 

Therefore,  Eq.  (XI  1)  can  only  be  valid  both  outside  the  singular 

layer  and  outside  the  island.  In  the  linear  regime,  where  §  is  infinitesimal, 

the  island  width  is  smaller  than  the  resistive  layer  width.  However  since 

L  is  small,  as  B  grows,  soon  Ax.  >  L  so  the  island  width  becomes  the 
c  is  c 
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m 


relevant  inner  layer  width.  Equation  (XI  1)  then  is  only  valid  for 
x  >  Axl8. 

We  now  show  how  the  quasl-llnear  evolution  of  a  single  tearing  mode 
may  be  carried  out.  If  the  variables  are  renormalized  in  dimensionless 
form  as 


p  ■  kx 


(a) 


x  -  fydt 


(b) 


the  equations  for  and 


sr 


(XI  3) 


(a) 


(b) 


(XI  4) 


*  §*  ?<gi>  +  yg/  /P<R1>\X  (O 

<V  V  V  ^  l  / 
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1 


(2  8/4<By>  p  -  0) 

>  3d  / 


To  get  the  second  form  of  Eq.  (XI  4c),  we  have  expanded  the  derivatives  and 

s2‘, 

used  Eq.  (XI  4a)  to  express  ^7 ■  .  Despite  the  apparent  complexity,  the 

second  form  in  Eq.  (XI  4c)  is  much  simpler  Lo  work  with  because  the  diffusion 

term  now  has  a  positive  diffusion  constant.  Also,  we  have  assumed  B^  is 

real.  Equation  (XI  4a)  is  the  equation  for  B  in  the  outer  region  (the  same 

X  '  *  1/2 
as  Eq.  VII  9)  except  that  now  |p|  >  where  p^s 

(The  additional  factor  of  ^7  between  this  and  Eq.  (VII  41)  3p 

comes  from  the  difference  between  sin  cos  and  exponential  notation). 

Equation  (XI  4b)  simply  says  that  the  fluctuation  at  p  *  p  grows  as  the 

IS 

linear  growth  rate.  Equation  (XI  4c)  is  the  quasi-linear  equation  for  (B  ) 

y 

which  is  valid  only  outside  the  island,  that  is  Ip  I  >  p 

is 

The  question  now  is  how  to  calculate  the  island  width,  which  depends 

upon  the  slope  of  < B^  at  x  -  0,  that  is, inside  the  singular  region.  Clearly 

it  is  necessary  to  make  some  assumption  concerning  the  current  and  field 

distribution  inside  the  singular  layer.  We  make  the  reasonable  assumption 

that  <By)  Is  linear  in  x  within  the  island.  Then  Eqs.  (4)  must  be  solved 

numerically  in  space  and  time,  coupled  with  a  calculation  of  the  island 

width.  The  procedure  is  as  follows. 

it  t  ■  0  start  with  a  field  profile,  which  we  will  specify  as  (  B  )  * 

B0  tanh  (p/La).  Here  Lg  is  the  scale  length  in  units  of  k~^.  To  advance 

from  t  to  x  +  At,  first  solve  Eq.  (XI  4a)  for  B  .  From  B  ,  calculate 

x  x 


A  = 


_i _  fHJ  _  SB,  I  ~)  -l 

s.ff/O  l  >r  Ip  if  l-p, ) 


SB, 

5* 


due  to  the  symmetry.  Note  that  now  4  is  calculated  across  the  island 
rather  than  as  a  discontinuity  at  x  -  0.  This  means  the  4  in  Eq.  (XI  5) 
is  proportional  to  the  field  energy  liberated  outside  the  island.  If 


(XI  5) 
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A  <  0,  the  plasma  has  become  stable  and  the  calculation  Is  finished  with 
the  <By)  *(  By(-r)).  If  A  >  0  the  plasma  is  still  unstable  so  that  By(p**pig,t) 
is  advanced  according  to  Eq.  (XI  4b)  and  (By)  is  advanced  according  to 
Eq.  (XI  c).  To  calculate  the  island  width  first  calculate  the  quantity 


r~r)/i  83(p, 


(XI  6) 


This  is  the  island  width,  determined  at  each  p,  assuming  the  field  inside 
the  island  is  linear  up  to  that  p  and  equal  to  (B^(p))  at  p.  The  actual 
island  width  is  then  determined  at  each  time  step  by  solving 


U  ( A  r*  '  f>  (xi  7 

and  the  solution  of  Eqs.  (XI  4a  and  c)  is  invoked  only  outside  the  island. 

The  field  inside  the  island  then  is  linear  in  p,  and  Eq.  (XI  7)  insures  that 

(  B^  is  continuous  at  p  =  Lig(p,t  +  At)« 

Equations  (XI  4a-c)  were  solved  numerically  by  Dr.  Barbara  Mellander 

using  this  assumption  for  the  field  inside  the  island.  The  initial  and  final 

magnetic  fields  for  four  choices  of  Lg  are  shown  in  Fig.  (XI  2a-d).  For 

weakly  unstable  plasma  (larger  L  ),  the  asymptotic  field  is  a  very  smooth 

s 

function  of  space.  However,  as  L  gets  smaller  there  is  a  larger  and  larger 

s 

discontinuity  in  derivative  at  the  island  edge.  Most  likely  this  means  that 
the  quasi-linear  theory  becomes  less  accurate  as  the  plasma  becomes  more 
unstable.  The  arrow  on  the  horizontal  axis  of  the  four  graphs  is  the  position 
of  the  maximum  of  B  at  t  ■  0.  For  weakly  unstable  plasmas,  the  final  island 
width  turns  out  to  be  smaller  than  this, for  strongly  unstable  plasmas,  larger. 
In  Fig.  (XI  3)  is  shown  both  the  final  island  width  and  final  island  width 
divided  by  L  as  a  function  of  L  .  In  Fig.  XI  4  is  shown  the  liberated 

8  S 

magnetic  energy  per  unit  area  and  the  liberated  magnetic  energy  per  unit  area 
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divided  by  as  a  function  of  scale  size.  To  summarize,  the  quasi- 

linear  evolution  of  the  outer  region  always  drives  a  tearing  mode  to  a  final 
stable  state  with  A  *  0  if  one  assumes  (B  )  is  linear  in  x  within  the  inner 

y 

region.  Also,  this  state  will  be  the  same  no  matter  what  is  the  dissipation 
mechanism  in  the  inner  region.  However,  the  quasi-linear  theory  almost 
certainly  gets  less  and  less  accurate  as  the  mode  gets  more  and  more  unstable. 

The  next  question  is  whether  there  are  nonlinear  effects  within  the 
island  itself  that  either  stabilize  the  mode  or  strongly  affect  its  behavior. 

This  problem  was  analyzed  by  Rutherford  (P.  Rutherford,  Phys.  Fluids  16,  1903  (1973) 
He  found  that  nonlinear  effects  in  the  island  did  not  stabilize  the  mode,  but 
did  significantly  reduce  its  growth  rate.  Here  we  give  a  very  much  simplified 
version  of  his  theory.  To  start,  recall  that  Eq.  (VII  16)  gives  the  power 


liberated  per  unit  area  in  the  outer  region  as 


y  a  «  g*  3* 

*3?  ” 


(XI  8) 


assuming  Bx  is  real.  This  power  is  dissipated  by  Ohmic  heating  within  the 
inner  region  of  width  L^.  That  is,  the  current  associated  with  the  dis¬ 
continuity  in  By  is  assumed  to  flow  essentially  uniformly  within  Lc,  so  that 

***  r  A 

this  current  is  given  by  Eq.  (VII  24),  J z  -  ^  Bx^x  *  °)  jp  Balancing 

c 

the  Ohmic  heating  with  this  power  liberated  yielded  the  approximate 


expression  for  the  growth  rate,  Eq.  (VII  25b). 

If  one  now  makes  the  reasonable  assumption  that  for  Ax 


..  •  i 


the  current  flows  throughout  the  island,  the  simple  analysis  leading  up 


to  Eq.  VII  25b,  is  replaced  with  ^ 


(XI  9) 
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or  for  long  time 


Thus  the  Island  growth  becomes  linear,  rather  than  exponential  In  time,  and 

the  growth  rate  goes  as  first  power  of  resistivity,  rather  than  some 

fractional  power.  Thus  once  Ax.  >  L  ,  the  island  growth  rate  does  not 

18  S 

stop,  but  it  does  drastically  slow  up,  even  if  A  does  not  change.  Hence 
there  are  two  complimentary  aspects  to  the  nonlinear  evolution  of  tearing 
modes.  First,  as  the  background  evolves,  A  is  reduced  and  eventually 
brought  to  zero  by  the  quasi-linear  evolution  of  the  outer  region  magnetic 
field  structure.  Second,  because  of  the  change  in  the  character  of  the 
inner  region  for  Ax^g  >  Lc>  the  growth  of  the  island  width  is  drastically 
reduced  and  it  proceeds  on  the  much  slower  resistive  diffusion  time  scale, 
the  same  time  scale  as  for  changes  in  equilibrium  quantities. 

Now  let  us  discuss  the  implication  of  this  for  minor  and  major 
disruptions  in  tokamak  plasmas.  As  was  discussed  in  Chapter  I,  a  minor 
disruption  is  a  sudden  burst  of  x-ray  activity  generally  associated  with 
a  measurable  disturbance  at  some  m  number,  usually  2  or  3.  The  time  scale 
for  this  x-ray  activity  in  PLT  was  typically  100  p  sec.  One  possible 
interpretation  is  that  tearing  modes  of  the  associated  helicity  are  excited. 
The  plasma  is  not  in  an  absolutely  steady  state  but  is  continuously  evolving 
due  to  for  instance  changes  in  the  external  circuit.  These  changes 
propagate  to  the  interior  of  the  plasma  in  a  resistive  skin  time,  about 
100  miliseconds  or  more.  As  the  interior  of  the  plasma  changes,  it  might 
evolve  to  an  unstable  state  with  A  >  0.  Then  a  tearing  mode  will  be 
excited.  Its  growth  time  will  be  very  small  compared  to  the  time  for  plasma 
evolution.  For  PLT,  Eq.  VII  25b  gives  a  growth  time  of  between  about  100  u 
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sec  and  1  milisec  depending  on  which  parameters  one  chooses.  The  theory 
developed  in  this  chapter  would  predict  that  the  plasma  then  evolves  to  a 
stable  state  at  which  either  A  =  0  for  complete  stability,  or  Ax^g  i 
for  drastic  reduction  in  growth  rate.  The  time  scale  would  be  several 
growth  times.  This  would  be  manifest  by  the  plasma  very  suddenly  changing 
its  equilibrium,  accompanied  by  an  x-ray  burst  and  the  onset  of  island 
structure. 

Now  let  us  discuss  whether  a  single  tearing  mode  could  lead  to  a 
major  disruption.  Imagine  that  at  some  radius  a  tearing  mode  is  excited. 

At  this  radius  an  island  or  chain  of  islands  is  formed.  While  the  tearing 
mode  evolves  quickly  compared  to  the  magnetic  diffusion  time,  it  evolves 
slowly  compared  to  the  Alfven  time.  Therefore  each  of  the  magnetic  surfaces 
within  the  island  is  an  MHD  equilibrium.  Hence  as  the  tearing  mode  evolves 
new  MHD  equilibria  form,  each  flux  surface  in  the  island  having  some  average 
of  the  pressure  of  the  two  or  more  original  flux  surfaces  which  reconnected 
to  form  it. 

If  the  island  grows  to  a  sufficient  size  that  it  touches  the  limiter, 
the  cold  plasma  at  the  limiter  will  flow  freely  along  the  field  and  end  up 
on  the  other  side  of  the  island,  deep  in  the  interior  of  the  plasma.  Not 
only  will  this  cold  plasma  get  to  the  interior,  but  in  all  likelihood  large 
numbers  of  metal  impurity  ions  from  the  limiter  will  also.  Radiation  from 
these  impurities  will  further  cool  the  plasma.  If  the  plasma  cannot  adjust 
to  this  sudden  cooling  of  the  interior,  it  may  be  that  a  major  disruption 
will  result.  Thus  it  is  likely  that  a  tearing  mode  can  produce  a  major 
disruption  in  a  tokamak  if  the  island  can  grow  until  it  reaches  the  limiter. 
This  hypothesis  seems  to  be  confirmed  by  the  results  on  Pulsator  discussed 
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in  the  introduction.  Here  an  island  is  artificially  induced  by  external 
coils,  and  when  the  island  touches  the  limiter,  a  major  disruption  did 
result. 

This  chapter  discussed  so  far  how  an  initially  tearing  mode  unstable 
plasma  in  slab  geometry  evolves  toward  a  final  state  having  a  periodic 
chain  of  islands.  We  will  close  the  chapter  with  a  discussion  of  the 
subsequent  evolution  of  this  island  chain.  If  the  original  equilibrium 
variation  was  in  the  x  direction,  this  time  asymptotic  periodic  chain  of 
islands  extends  in  the  y(horizontal) direction.  Each  island,  of  course, 
represents  a  current  filament  in  the  z  direction.  Since  these  currents 
are  all  in  the  same  direction,  they  attract  each  other.  The  plasma  is  in 
equilibrium  because  an  island  feels  equal  and  opposite  attractive  forces 
from  the  island  to  the  right  and  the  island  to  the  left. 

Now  consider  what  happens  if  the  island  is  given  a  rightward  displace¬ 
ment.  First  of  all  the  attractive  force  of  the  two  nearby  islands  is  greater 
than  the  attractive  force  to  their  now  more  distant  neighbors.  Thus  the 
nearby  islands  tend  to  attract  each  other  and  ultimately  coalesce  so  if  this 
effect  is  dominant,  the  island  chain  is  unstable.  This  attractive  force, 
however,  may  be  balanced  by  a  repulsive  force.  As  the  two  islands  move 
toward  each  other,  the  flux,  which  is  frozen  into  the  flow,  must  be  compressed 
between  the  two  islands.  Thus  the  magnetic  pressure  increases  there  and 
this  forces  the  two  islands  apart.  If  this  latter  effect  dominates,  the 
plasma  is  stable,  at  least  to  this  type  of  displacement.  The  calculation 
of  the  stability  of  this  plasma  is  very  complicated  due  to  the  two  dimensional 
nature  of  the  equilibrium.  One  calculation  has  been  made  by  J.  Finn  and 
P.  Kaw  (Phys.  Fluids  20,  72  (1977))  and  they  found  that  the  attractive  force 
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generally  dominates  so  that  the  island  chain  is  unstable  to  coalescence 
of  the  individual  islands. 

A  numerical  simulation  of  the  evolution  of  this  equilibrium  was  done 
by  P.  Prichett  and  C.  Wu  (Phys.  Fluids  22,  2140  (1979)).  At  t  =  0  an  MHD 
equilibrium  of  a  periodic  chain  of  islands  was  set  up,  the  flux  surfaces 
of  which  are  shown  in  Fig.  (XI  5).  The  simulation  had  periodic  boundary 
conditions  in  the  horizontal  (y)  direction,  each  periodicity  length 
containing  two  islands.  If  the  plasma  had  zero  resistivity,  the  islands 
displaced  toward  each  other,  but  ultimately  stopped  when  the  magnetic 
field  compressed  between  them  became  large  enough  to  repel  their  coalescing 
motion.  The  flux  surfaces  at  two  later  times  are  shown  in  Fig.  (XI  6). 
However,  if  resistivity  is  present,  the  field  structure  can  change  its 
topology  and  the  islands  can  coalesce.  The  flux  surfaces  for  a  simulation 
with  n  ^  0  are  shown  in  Fig.  (XI  7).  Clearly  the  magnetic  reconnection 
proceeds  to  completion  and  the  two  initial  islands  merge  to  form  a  single 
island. 

Thus  the  nonlinear  evolution  of  a  tearing  mode  unstable  plasma  can 
be  quite  complex.  There  are  at  least  three  processes  which  can  play  a 
role,  quasi-linear  evolution  of  the  plasma  in  the  outer  region,  reduction 
in  growth  due  to  the  dynamics  of  the  plasma  in  the  island,  and  ultimately 
coalescence  of  the  resulting  island  chain. 
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XI  5)  The  Initial  Island  structure  examined  by  Prlchett  and  Wu 


XII.  Island  Overlap  and  the  Onset  of  Stochaatlcity 

In  the  last  chapter  we  discussed  the  response  Oi.  the  plasma  to  a 
single  tearing  mode  In  slab  geometry.  Although  the  magnetic  topology  Is 
changed  by  a  tearing  mode,  magnetic  surfaces  still  do  exist.  In  this 
section  ve  Investigate  under  what  circumstances  unstable  MHD  fluctuations 
can  destroy  magnetic  surfaces.  For  simplicity  we  consider  cylindrical 
geometry  to  lowest  order  in  tokamak  ordering,  so  that  B  is  two  dimensional 
in  r  and  9  and  is  constant  and  very  large. 

Hence  B  follows  from 
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(XII  1) 


where  A  is  the  z  component  of  the  vector  potential.  The  equation  for  the 
field  line  is 

Ar  _  Br  _  JL 

<sl *■  j  r  &Z.  (XII  2) 

There  are  at  least  two  possible  ways  to  display  the  solution  of  Eqs.  (XII  2). 
First  of  all,  the  solutions  could  be  projected  to  z  ■  0,  the  result  being 
a  curve  in  the  SB  plane.  Secondly,  if  the  system)  is  periodic  in  z  with  period 
2dR  (obviously  like  a  tokamak  with  radius  R),  the  intersection  of  the  field 
line  with  the  planes  z  ■  2imR  could  be  displayed  on  the  z  ■  0  plane.  The 
result  now  is  not  a  curve,  but  a  series  of  points  which  may  either 
lie  on  a  curve,  or  else  fill  an  area. 
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Making  use  of  Eq.  (XII  1),  it  is  clear  that  Eqs.  (XII  2)  are 
Hamiltonian  in  form,  with  z  playing  the  role  of  t,  and  A  the  role  of  the 
Hamiltonian.  Since  A  depends  on  z,  the  Hamiltonian  is  not  a  constant  of 
the  motion.  However  since  the  motion  of  the  field  line  does  follow  from 
a  Hamiltonian , the  transformation  from  r(z  ■  0),  0  (z  *  0)  to  r(z),  9  (z)  is 
an  area  preserving  transformation.  There  has  been  a  tremendous  amount  of 
recent  work  on  whether  such  area  preserving  transformations  have  an  ordered 
or  stochastic  nature.  (A  good  introduction  and  review  can  be  found  in 
J.  Ford,  Fundamental  Problems  in  Statistical  Mechanics  III,  E.  Cohen  ed 
1974).  The  basic  motivation  behind  this  work  is  to  see  when  statistical 
mechanics  is  valid  for  an  Isolated  Hamiltonian  system.  This  work  can  be 
quite  mathematically  abstract.  For  instance  it  has  only  recently  been 
proven  rigorously  that  a  hard  sphere  gas  obeys  the  ergodic  theorm.  To 
quote  Ford,  this  proof  occupies  "about  one  hundred  journal  pages  and  will 
involve  such  concepts  as  discontinuous  transverse  foliations,  completely 
positive  Kolmogorov  entropy,  and  the  like."  Of  course  it  is  far  beyond 
the  scope  of  this  work  to  more  than  scratch  the  surface  of  the  mathematical 
theory  of  ergodic  behavior.  Principally  we  are  Interested  in  whether  the 
intersection  of  field  lines  with  z  ■  2imR  (or  what  we  will  call  the  equivalent 
z  ■  0  plane)  fill  this  plane,  or  some  region  of  it  ergodically,  or  whether 
the  Intersections  like  on  a  well  defined  curve. 

One  thing  is  very  easy  to  show.  Namely  if  there  does  exist  a  constant 
of  motion  say  x(rfi  ,z)  -  constant,  then  the  intersections  of  field 

lines  with  the  equivalent  z  -  0  plane  line  on  a  well  defined  curve;  in 
other  words  the  motion  of  the  field  line  is  not  ergodic.  To  show  this. 
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first  note  that  since  the  system  is  periodic  in  z, 

V  (r,0,z)*/  (C 

Therefore  on  the  equivalent  z  *  0  plane, 

/  (r, 9,  z=°)  = 


(XII  3) 


(XII  4) 


so  that  the  intersections  of  a  field  line  with  this  plane  lie  on  the 
curve  defined  by  Eq.  (XII  4). 

We  will  now  prove  that  if  the  perturbed  magnetic  field  has  helical 
symmetry,  that  is  the  dependence  on  6  and  z  is  through  the  combination 
t  ■  me+  kz,  then  a  constant  of  motion  exists.  To  do  so,  note  that  V*B  *  0, 


2-r  3r  +  ■»>  —  +  r*  — 

?rrU)r  ?r  ir 

Therefore  B  can  be  expressed  in  terms  of  a 


=  O 

'helical  flux  function' 


(XII  5) 
as 
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The  condition  that  Eqs.  (XII  6)  have  a  solution  for  if  B  is  given,  is 

simply  that  the  divergence  of  B  is  equal  to  zero.  Notice  that  in  general, 

Eq.  (XII  6)  cannot  be  used  to  solve  for  B  in  terms  of  ^  since  it  consists 

of  only  two  equations  for  the  three  components  of  B.  (If  tokamak  ordering, 

or  B  ■  constant,  is  valid,  then  Eqs.  (XII  a  and  b)  also  solve  for  B  if 
z 

is  given.)  However  the  information  in  Eq.  (XII  6)  is  sufficient  to  show 
that  a  constant  of  motion  exists.  It  follows  directly  from  Eqs.  (XII  6) 
that 

8*  V  —o  '(xii  7) 

so  that  B  lies  in  the  surfaces  of  constant  In  other  words  ^  is  constant 
along  the  line  of  force. 
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Now  consider  a  magnetic  field  with  a  helical  perturbation  specified  by 


[j)  (r,t)  ~  +  €(r)  CoQ'  ^ 


(XII  8) 


and  for  simplicity,  take  also  k 


Note  that  the  derivation  of  the  Oi  term 

o 


vanishes  at  radial  position  r  given  by 
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vrt 


(XII  9) 


Near  this  singular  surface  the  surfaces  of  constant  are  given  by 


r-r. 


/  K-  cos  r 

"  l  *  *7™  J 


i  K  y-  V'/CO 

so  that  the  surfaces  of  constant  i|*  look  as  shown  in  Fig.  (XII  1). 
may-imnm  radial  island  (ie  when  K  =  e)  with,  from  center  to  boundary, 

,  (  e<o  Vi 

4r.-s=  i\yr) 


(XII  10) 
The 

is  given  by 


(XII  11) 


returning  to  standard  cylindrical  geometry,  the  surfaces  of  constant  \p  are 
shown  in  Fig.  (XII  2)  for  the  case  of  s  «  3,  The  projection  onto  the 
equivalent  z  -  0  plane  is  shown  in  Fig. (XII  3).  Thus  for  the  case  of  a 
single  perturbation,  the  helical  symmetry  allows  us  to  prove  the 
existence  of  constant  of  motion*  Therefore  the  field  line  is  not  ergodic. 


The  obvious  question  next  is  what  happens  if  the  helical  symmetry  is 
broken,  for  instance  by  the  presence  of  an  additional  perturbation  with  k 
still  equal  to  -'Vfe  but  a  different  asimuthal  wave  number  m  •  Now  there  is 
no  single  helical  symmetry,  and  apparently  no  constant  of  motion.  This 
problem  has  been  investigated  extensively  by  numerical  simulation.  While 
a  fact  of  ergodicity  still  eludes  a  rigorous  mathematical  proof,  the  results 
of  many  numerical  simulations  and  analytic  theories  does  point  toward  a  very 
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reasonable  hypothesis.  Namely,  if  the  island  widths  of  the  m  and  m 

perturbation  are  small  compared  to  the  inter-island  separation,  well 

/ 

define  islands  exist  around  the  q  =  m  and  q  =  m  points.  However  the 

separatrices  may  not  be  distinct  lines  but  rather  thin  ergodic  regions. 

! 

Between  the  q  =  m  and  q  =  m  points,  the  field  line  projections  are 
mostly  as  in  Fig.  (XII  3)  except  there  may  be  small  regions  of  ergodicity 
and  small  chains  of  secondary  islands,  as  indicated  in  Fig.  (XII  4a)  for 
an  m  =  2  and  m  =  3  perturbation  with  all  islands  having  the  same  width. 

As  the  island  width  increases,  the  ergodic  regions  increase  around  the 
separatrix  as  shown  in  Fig.  (XII  4b).  Finally,  when  the  island  edges 
overlap,  most  of  the  region  between  the  q  =  2  and  q  =  3  surface  is  ergodic, 
except  for  perhaps  a  small  region  in  the  center  of  each  island,  as  shown 
in  Fig.  (XII  4c).  Thus  the  condition  of  island  overlap  is  almost  universally 
accepted  to  be  the  condition  for  ergodic  behavior  between  the  relevant 
rational  surfaces. 

Let  us  see  what  this  implies  for  the  case  of  a  m  =  2  and  m  =  3 
perturbation.  The  distance  between  the  two  rational  surfaces  q  *  2  and 
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Thus  in  a  system  with  strong  shear  (large  q ") ,  island  overlap  is  accom¬ 
plished  with  smaller  perturbed  field.  Each  island  width  is  smaller, 
but  the  distance  between  neighboring  rational  surfaces  is  also  less. 

Consider  now  the  implication  of  this  for  a  tokamak  plasma.  Imagine 
that  the  plasma  is  initially  unstable  to  both  an  m  =  2  and  m  =  3  tearing 
mode.  In  the  early  stages,  when  the  island  width  is  small  compared  to 
island  separation,  the  equilibrium  must  readjust  to  the  presence  of  the 
fluctuations,  but  since  magnetic  surfaces  exist,  equilibrium  also  at  least 
exists.  However  when  the  perturbations  grow  to  such  an  amplitude  that 
these  islands  overlap  MHD  equilibrium  is  suddenly  (ie  within  a  growth 
time)  lost  in  a  very  large  region  of  the  plasma,  and  the  only  possible 
pressure  profile  is  p  =  constant  in  this  extensive  ergodic  region.  If 
this  sudden  and  violent  re-arrangement  of  the  pressure  profile  is 
sufficient  to  cause  a  major  disruption  of  a  tokamak  discharge,  then 
island  overlap  is  one  possible  cause  of  such  major  disruptions. 


XII  3)  The  projection  of  the  island  structure  onto  the 
equivalent  Z  ■  0  plane  for  m  *  3 


XII  4)  The  projection  of  the  island  structure  onto  the  equivalent  Z  ■  0 
plane  for  the  case  of  equal  m  *  2  and  m  ■  3  perturbations.  The  dotted 
regions  represent  the  region  of  stochastic  field  lines.  As  the  amplitude 
of  the  perturbations  increase,  more  of  the  region  is  stochastic. 


XIII.  The  Taylor-Woltjer  Theory  of  Spontaneous  Field  Reversal  and  Current  Limitation 


The  past  three  chapters  have  discussed  nonlinear  theory  with  reference  to 
some  particular  instability.  Other  authors,  however,  have  looked  at  the  non¬ 
linear  motion  of  an  unstable  plasmas,  but  without  specific  reference  to  a 
particular  instability.  In  this  chapter  and  the  next,  we  survey  some  of  this 
work.  This  chapter  discusses  the  theory  of  Woltjer  (Proc.  Natl.  Acad.  Sci. 

U.S.  44,  489  (1958))and  Taylor  (Phys.  Rev.  Lett.  33,  1139  (1974)),  for  spontaneous 
reversal  of  the  axial  field  and  of  current  limitation  in  a  reversed  field  pinch. 
The  fundamental  assumption  in  this  work  is  that  a  plasma  will  relax  to  its 
lowest  energy  state  consistent  with  all  the  constraints  on  it.  Taylor's  main 
initial  assumption  is  that  the  thermal  and  flow  energy  of  the  plasma  is  much 
less  than  the  magnetic  energy,  /d3r  B2/8tt.  This  (unfortunately)  is  true  in 
current  tokamaks  and  pinches,  although  the  ultimate  hope  is  certainly  to  attain 
a  high  beta  plasma.  The  fact  that  the  pressure  is  vanishingly  small,  of  course, 
imposes  severe  constraints  on  the  MHD  equilibrium.  Consider  for  a  moment  the 
MHD  equilibrium  of  a  high  aspect  ratio  (^/a  »  1)  torus;  that  is  the  topology 
is  that  of  a  torus,  but  all  vector  equations  can  be  expressed  in  cylindrical 
co-ordinates.  This  equilibrium  now  becomes 


£(v  *?)*?=■  Vp-o 


(XIII  1) 

Equation  (XIII  1)  simply  says  that  V.  x  B  is  everywhere  parallel  to  B,  or 


Vxs=  xnr)I 

"  (XIII  2) 

where  A  is  some  scalar  function  of  position.  Taking  the  divergence  of 

both  sides  of  Eq.  (XIII  2)  and  making  use  of  the  fact  that  V*B  ■  0,  we  find 
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so  that  A  is  constant  along  every  flux  tube,  flux  surface,  or  flux  volume; 
whichever  is  appropriate. 


Taylors  theory  then  gives  a  prescription  for  finding  A  and  also 
determining  the  resulting  boundary  conditions  for  Eq.  (XIII  2).  It  does 
this  by  minimizing  the  magnetic  energy  subject  to  appropriate  constraints. 
The  first  order  of  business  then  is  to  find  these  constraints;  minimizing 
magnetic  energy  subject  to  no  constraints  gives  only  the  trivial  solution 
B  -  0. 

The  configuration  is  a  toroidal  pinch  bounded  by  a  conductor.  However, 
as  shown  in  Fig.  (XIII  1),  there  is  a  poloidal  slit  with  a  voltage  source 


across  it.  At  t  *  0,  J  ■  0  and  there  is  only  a  toroidal  field  specified  by 


8*  8,-  £, 
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Then  at  t  ■  0,  a  voltage  is  pulsed  across  the  gap  for  a  time  6 t.  For 
later  times,  the  gap  is  shorted  out  (crowbarred).  The  plasma  then  responds 
to  this  Induced  voltage,  and  relaxes,  presumably  due  to  instabilities,  to 
some  final  minimum  energy  state.  Since  the  plasma  is  surrounded  by  a 
perfect  conductor,  no  toroidal  flux  can  get  in  or  out,  so  the  toroidal 
flux  is  conserved  and  has  value  V  ■  irB^a2.  However  because  of  the  poloidal 
slit,  poloidal  flux  can  go  in  or  out  during  the  time  6t,  so  it  is  not 
conserved.  Thus  the  first  constraint  is  that  ¥  is  constant.  As  we  will 
see,  this  determines  the  boundary  conditions. 

A  second  constraint  follows  from  the  fact  that  the  magnetic  flux  is 


frozen  into  the  fluid  flow.  We  will  now  show  that  this  Implies  that  as 
the  fluid  relaxes,  K  ■ ^A»Bd3r  is  constant  after  the  voltage  pulse  is  shut 
off.  Since  B  ■  7  x  A,  Maxwell's  equation  for  (Eq.  II  8)  can  be 


integrated  once  to  give 

2i  .  = 


(XIII  5) 
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where  w  is  any  scalar  function.  The  presence  of  *  on  the  right  hand  side 
of  Eq.  (XIII  5)  of  course  accounts  for  the  guage.  Taking  the  dot  product 
of  Eq.  (XIII  5)  with  B  and  integrating  over  space  gives  the  result 


f/A8-0-  JJfi'rB- If: 
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The  surface  integral  arises  because  V*B  ■  0;  it  vanishes  because  B  is 


everywhere  normal  to  the  bounding  surface.  Thus 

iff/S'S-l-  #'• 
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Nov  examine  the  two  terms  on  the  extreme  right  hand  side  of  Eq.  (XIII  7). 

In  a  perfectly  conducting  plasma  E  is  perpendicular  to  B,  so  the  first 
term  vanishes.  To  evaluate  the  second,  we  need  the  components  of  E  and  A 
only  in  the  plane  of  the  bounding  surface.  After  the  system  is  crowbarred, 
E  is  zero  in  the  plane  of  the  surface  (since  the  surface  is  a  perfect 
conductor),  so  A»B  is  constant. 


For  the  time  6t,  during  which  the  system  is  pulsed,  E  *  V  /<5  is  the 

z  o 


width  of  the  slit,  so  that 
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Thus  the  rate  of  change  of  K  is  simply  proportional  to  the  product  of 
voltage  times  toroidal  flux.  Since  K  (t  *  0)  ■  0,  as  is  obvious  from 
Eq.  (XIII  4),  then 


K- 
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so  that  after  time  fit,  K  is  proportional  to  the  Volt-seconds  stored  in  the 
external  circuit. 


The  problem  now  is  to  determine  the  state  that  the  plasma  relaxes  to 
after  the  voltage  source  is  turned  off.  The  basic  assumption  is  that  the 
plasma  relaxes  (presumably  by  instability)  to  a  state  which  minimizes  the 
magnetic  energy,  subject  to  the  constraint  that  K  is  constant.  Using  the 
method  of  Lagrange  multipliers,  this  means  minimizing 
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where  p  is  the  Lagrange  multiplier.  It  is  a  simple  matter  to  show  that 
the  right  hand  term  in  Eq.  (XIII  10)  is  guage  invariant  by  letting 
A  -*■  A  +  Vy»,  integrating  by  parts,  and  making  use  of  the  fact  that  B  Is 
parallel  to  the  surrounding  conducting  wall.  To  minimize  the  expression 
in  Eq.  (XIII  10),  let  B  -*•  B  +  fiB,  A  -*■  A  +  6A  (of  course  6B  *  V  x  6A),  and 
set  terms  linear  in  6A  and  6B  equal  to  zero,  so  that 


Jfjj\  [sB-SB  -m  +  A-JB)]  =  0 


(XIII  11) 


Writing  B  and  6B  in  terms  of  A  and  5A  and  performing  various  partial 
integrations,  it  is  not  difficult  to  show  that  Eq.  (XIII  11)  reduces  to 


(XIII  12) 


where  the  second  integral  in  Eq.  (XIII  12)  is  over  the  bounding  conducting 
surface.  The  next  step  is  to  prove  that  this  surface  integral  vanishes. 

To  start,  consider  the  integral^A'di  around  any  closed  curve  on  the 
surface.  This  integral  is  just  the  magnetic  flux  linking  the  closed  curve. 


There  are  two  possibilities.  If  the  curve  links  the  torus  in  the  poloidal 
direction,  the  integral  is  Y,  the  toroidal  flux.  Otherwise  the  integral 
vanishes.  Since  all  variations  in  B  leave  the  flux  invariant. 


(XIII  13) 


Therefore  on  the  bounding  surface  6A  -  Vgg  where  is  a  two  dimensional 
gradient  within  the  surface,  and  g  is  a  scalar  function  defined  on  the 
surface.  Since  6A  normal  to  the  surface  gives  no  contribution  to  the 
second  term  in  Eq.  (XIII  12),  it  can  be  written  as 


(XIII  14) 
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The  second  term  on  the  right  of  Eq. 
over  a  closed  surface  of  the  curl 
can  be  expressed  as 


(XIII  14)  vanishes  since  the  integral 
of  any  vector  vanishes.  The  first  term 


(XIII  15) 


Because  the  normal  component  of  both  J  and  B  are  equal  to  zero  on  the 
conducting  surface,  the  right  hand  side  of  Eq.  (XIII  15)  is  zero.  Hence 
setting  the  constrained  variation  of  energy  equal  to  zero,  Eq.  (XIII  12) 
reduces  to  the  simple  result 


Vx6 
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Comparing  Eq.  (XIII  16)  to  Eq.  (XIII  2)  we  see  that  the  only 
difference  is  that  p  is  a  constant,  whereas  X  is  constant  over  a  flux 
tube  or  surface,  but  otherwise  can  vary  in  space.  Thus  minimizing  the 
energy  subject  to  the  constraint  K  ■  constant,  picks  out  a  particular 
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MHD  equilibrium  out  of  all  those  which  are  possible  according  to  Eq.  (XIII  2). 
Let  us  now  examine  how  this  selection  comes  about. 

In  deriving  Eq.  (XIII  16),  the  only  boundary  conditions  used  were  that 
the  components  of  B  and  J  normal  to  the  conducting  surface  both  vanish.  In 
ideal  MHO,  this  is  not  only  true  at  the  conducting  wall,  but  is  also  true  of 
each  flux  surface,  as  shown  in  Chapter  II.  Therefore,  one  could  equally  well 
minimize  the  magnetic  energy  between  two  flux  surfaces,  subject  to  the 
constraint  that  K  is  constant,  between  them,  and  derive  Eq.  (XIII  2)  in 
exactly  the  same  way  as  Eq.  (XIII  16)  was  derived.  Imagine  that  at  t  *  0 
the  flux  surfaces  are  circular,  as  shown  in  Fig.  (XIII  2a),  but  V  x  B  j*  A(Y)B 
(Y  denotes  the  flux).  As  the  plasma  releases  its  excess  magnetic  energy  and 
violently  evolves  toward  a  final  state  described  by  Eq.  (XIII  2),  these  flux 
surfaces  will  naturally  become  very  contorted,  but  must  maintain  their 
topological  structure  because  the  flux  is  frozen  into  the  flow,  as  discussed 
in  Chapter  II.  For  instance  at  t  *  00  infinity,  these  three  flux  surfaces 
might  evolve  toward  that  shown  in  Fig.  (XIII  2b).  In  this  final  state 
V  x  B  *  A(Y)B  on  each  flux  surface,  but  A  can  vary  from  one  flux  surface  to 
another.  Notice  that  there  are  places  in  Fig.  (XIII  2b)  where  flux  surfaces 
are  forced  together.  For  instance  at  the  star,  all  three  flux  surfaces  are 
forced  together.  Of  course  all  flux  surfaces  initially  between  1  and  3  are 
forced  together  at  this  point  also. 

Now  imagine  that  there  is  a  very  small  plasma  resistivity,  so  that  the 
field  lines  can  reconnect.  For  instance  if  resistivity  were  suddenly  turned 
on  for  the  configuration  in  Fig.  (XIII  2b),  surely  all  flux  surfaces  between 
1  and  3  would  reconnect  with  one  another  near  the  star,  and  other  reconnection 
would  occur  at  other  places  also.  Of  course  reconnection  changes  topology 
of  B,  but  only  has  a  small  effect  on  the  magnitude  of  B  in  most  places. 
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Therefore,  If  this  reconnection  Is  allowed,  it  makes  no  sense  to  talk 


about  magnetic  energy  or  constraints  between  two  flux  surfaces!  However 
all  relevant  boundary  conditions  can  still  be  applied  on  the  conducting 
wall,  so  that  minimizing  the  total  plasma  energy,  subject  to  the  constraint 
that  K  is  constant  is  still  valid  for  the  entire  plasma  even  if  not  for  each 
initial  flux  surface.  Thus  inherent  in  the  reduction  of  Eq.  (XIII  2)  (which 
could  have  also  been  derived  by  minimizing  the  energy  on  each  flux  surface) 
to  Eq.  (XIII  16)  is  the  assumption  that  in  the  violent  relaxation  of  the 
plasma  to  its  final  state,  a  small  amount  of  resistivity  is  allowed,  so 
that  topology  of  the  field  is  destroyed. 

Now  let  us  find  the  solution  to  Eq.  (XIII  16).  If  B  has  cylindrical 
symmetry,  it  is  a  simple  matter  to  show  that 

8,  --  8.  r.  (m  r) 
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Be  = 


is  a  solution  of  Eq.  (XIII  17),  and 
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(XIII  18) 
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Here  a  is  the  minor  radius  and  A has  an  appropriate  constant  added  to  it 
so  that  the  poloidal  magnetic  flux  through  the  hole  in  the  torus  vanishes. 
The  next  problem  is  to  relate  p  and  Bq  to  physical  parameters.  The  axial 
magnetic  flux  ¥  is  just  the  line  integral  of  the  poloidal  vector  potential, 
so 
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The  simplest  other  relation  is  that  between  p  and  the  pinch  ratio  c 

where  B^  is  the  initial  axial  bias  field  at  time  t  ■  0.  Relating 
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to  the  flux  and  the  current  to  ^  (r  ■  a),  we  find 


«B;C 

where  I  Is  the  current  in  CGS  units.  The  quantity  K  Is  related  to  the 


M, 

A 


stored  Volt  seconds  by  Eq.  (XIII  9).  Using  the  expression  for  the  fields 


given  In  Eq.  (XIII  17  and  18),  one  can  show  that 

K  K  f  mo.(s'(u*)  +  J 

o-L 


(XIII  21) 


3T  a  L  0 

which  defines  p  in  terms  of  initial  flux  and  stored  Volt  seconds. 

The  remarkable  thing  about  Eq.  (XIII  17)  is  that  it  predicts  that 
the  toroidal  field  reverses  direction  whenever  pa  >  2.4,  the  position  of 


the  first  zero  of  J  ,  or  for  pinch  parameter 


>  l.-z 
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a.  B;  C 

This  is  in  reasonable  agreement  with  the  results  of  pinch  experiments  as 
shown  in  Eq.  (I  ). 

So  far,  we  have  only  considered  solutions  of  Eq.  (XIII  16)  which  have 
cylindrical  symmetry.  These  are  also  solutions  which  are  not  symmetric. 
Expressing  the  azimuthal  and  axial  dependence  of  the  non-symmetrlc  solution 
to  Eq.  (XIII  16)  as  exp  i(kz  +  n6  ),  and  using  7»B  ■  0,  it  is  not  difficult 


to  show  that  Bz  satisfies  Bessels  equation,  so  that  the  solutions  are 


-  I 
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with  y2  =  (p2  -  k2)r2.  Of  course  Eq.  (XIII  23)  is  only  a  valid  solution 


to  the  equation  if 


(XIII  24) 


This  then  imposes  a  relation  between  p  and  k  which  must  be  satisfied. 

Since  Eq.  (XIII  16)  is  linear,  the  magnetic  field  can  be  the  symmetric 
state,  Eq.  (XIII  17),  plus  any  linear  combination  of  individual  solutions 
in  Eq.  (XIII  23).  The  condition  this  solution  must  satisfy  is  that  the 
total  K  and  toroidal  flux  is  given  and  that  B^(r  =  a)  =  0.  Since  the 
solutions  given  in  Eq.  (XIII  23)  have  zero  toroidal  flux,  this  toroidal 
flux  is  determined  entirely  by  the  symmetric  state.  The  quantity  K  then 
is  a  summation  over  contributions  from  each  solution  in  Eq.  (XIII  23)  plus 
the  contribution  from  the  symmetric  state. 

The  question  now  is  which  state  the  plasma  picks  out.  Of  course  the 
state  is  that  having  minimum  magnetic  energy.  Imagine  now  that  there  are 
two  solutions  to  Eq.  (XIII  16)  having  two  different  values  of  p  and  which 
satisfy  all  of  the  appropriate  boundary  conditions.  We  now  show  that  the 
minimum  energy  state  is  the  state  having  minimum  p.  The  energy  is  given  by 


W  --  JJpV  B-B  :  J/A  B-VWt 
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Writing  B  in  terms  of  A  in  Eq.  (XIII  16)  and  integrating  once,  it  is  easy 
to  show  that 


V*  /  -  ^  +  v  ^ 
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where  is  any  scalar  function.  Thus 

W=  jjJ®-  A^r '  mK 


(XIII  27) 


The  v  term  vanishes,  as  one  can  show  by  integrating  it  by  parts.  Therefore, 
since  K  is  constant,  the  minimum  energy  state  is  the  state  with  minimum  p. 

One  possible  state  for  the  plasma  is  simply  the  symmetric  state  of 
Eq.  (XIII  17).  However  if  a  helically  perturbed  state  with  smaller  p  can 
be  found,  the  plasma  should  relax  to  it.  Taylor  (J.  B.  Taylor,  Plasma 
Physics  and  Controlled  Fusion  Research,  1975,  IAEA  Vienna)  has  performed 
a  detailed  investigation  and  found  that  the  smallest  p  for  which 
Eq.  (XIII  24)  can  be  satisfied  is  pa  =  3.11  for  a  helical  perturbation 
having  m  =  1  and  ka  =  1.25.  Since  all  other  helical  states  have  higher  p, 
and  higher  energy,  we  need  only  consider  this  one. 

Now  review  how  the  plasma  state  is  set  up.  At  t  =  0,  there  is  some 
toroidal  flux  and  the  system  is  pulsed  with  a  certain  number  of  volt-seconds 
(ie  K).  A  current  (proportional  to  p)  is  thereby  induced  in  the  plasma. 

The  natural  expectation  is  that  I  (or  equivalently  p),  is  a  monotonically 
increasing  function  of  K  (the  volt  seconds).  However  for  a  sufficiently 
large  K,  p  will  be  equal  to  3.11.  As  K  is  further  increased  the  current 
will  no  longer  increase.  The  reason  is  that  a  lower  energy  state  exists 
which  is  a  linear  combination  of  the  cylindrically  symmetric  state  and 
helical  state  with  p  =  3.11.  The  toroidal  flux  is  determined  by  the 
cylindrical  state,  and  the  relative  amplitudes  of  the  cylindrical  and  helical 
state  are  specified  by  K.  Thus  after  p  *  3.11,  increasing  the  Volt  seconds 
does  not  increase  the  current,  but  rather  increases  the  amplitude  of  the 


helical  displacement.  The  current  limitation  predicted  here, 

a  c 


(XIII  28) 


is  also  in  reasonable  agreement  with  what  is  measured  in  pinch  experiments 
as  given  in  Eq.  (I  ).  Thus  this  relatively  simple  theory  gives  quite 


good  agreement  with  experiment  on  the  two  really  amazing  features  of 
reverse  field  pinches,  spontaneous  field  reversal  and  current  limitation. 

To  conclude,  it  is  worth  noting  that  Gibson  and  Whiteman  (Plasma 
Phys.  10,  1101  (1968))have  examined  the  stability  of  the  cylindrically 
symmetric  state  to  tearing  modes.  They  have  shown  that  this  transition 
from  a  cylindrically  symmetric  state  to  a  helically  perturbed  state 
corresponds  to  the  marginally  stable  point.  Of  course  this  is  to  be 
expected,  since  as  was  shown  in  Chapter  VII,  the  tearing  mode  first  goes 
unstable  when  the  magnetic  structure  has  a  neighboring  equilibrium  at 
the  same  energy.  The  theory  in  this  chapter  however  is  quite  different 
from  linear  stability  theory  in  that  growth  rates  are  not  calculated, 
but  the  amplitude  of  the  helical  perturbation  is. 
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XIV.  Kadomtsev's  Theory  of  Internal  Disruptions  and  Introduction  to  Numerical 
Simulations 


This  final  chapter  begins  with  Kadomtsev's  theory  for  m  =  1  internal 
disruptions  in  a  tokamak*  Amazingly  enough,  he  is  able  to  develop  this 
theory  without  ever  getting  involved  in  the  complications  of  the  m  =  n  =  1 
internal  kink-tearing  mode.  Preliminary  to  discussing  this  theory ,it  is 
necessary  to  review  the  toroidal  flux  t|t first  introduced  in  Chapter  XII.  Since 
all  dependence  on  Ooccurs  in  the  combination  r  =  +  kz,  it  was  shown  that 

the  fields  could  be  derived  from  Eq.  (XII  6).  In  lowest  order  tokamak 
ordering  where  B^  =  constant,  this  is  sufficient  to  solve  for  and  B^ 

consider  only  lowest  order  tokamak  ordering  here).  It  was  also  shown  that  ^ 
actually  is  a  flux,  it  is  simplest  to  relate  it  to  the  vector  potential.  By 
making  use  of  the  fact  that  B  =  Vx  A,  it  is  not  difficult  to  show  that  Eqs. 


(XII  6)  follow  if 

Ip* 


(XIV  1) 


To  continue,  we  will  relate  .Jf to  the  flux  through  a  helical  ribbon  at  radius 
r  and  pitch  defined  t  =  kz  +  m  0  =  constant.  The  flux  through  it  is  of 
course  just  the  line  integral  of  A  around  it.  A  unit  vector  parallel  to 
the  edge  of  the  ribbon  is 


making  use  of  lowest  order  tokamak  ordering.  Therefore  the  magnetic  flux 
through  the  ribbon  is  given  by  2ir  R  \f>,  where  R  is,  as  usual,  the  major 
radius  of  the  torus.  Hence  ip  is  indeed  the  helical  flux. 
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m  =  1  and  examine 


We  now  further  specialize  to  the  case  of  k  ® 

the  form  of  <|>  (r).  As  shown  in  Chapter  XII,  ip'  =  0  at  the  resonant  surface 
rB 

2 

where  q  =  =  1.  The  second  derivative  of  ^  at  the  resonant  surface  is 

a/*t/  _  ol  __  _0£  __  _  g  <  Q 

Jr  Ax  R  Jr  (xiv  2) 

where  we  have  assumed  the  normal  tokamak  profile  having  q  >  0.  Thus  ip 
has  a  maximum  at  resonant  surface  and  has  the  r  dependence  shown  in  Fig. 

(XIV  la). 

Let  us  now  examine  how  magnetic  reconnection  could  proceed  in  an  m  =  1 
internal  disruption.  The  magnetic  flux  through  a  ribbon  of  width  dr  is 
^7  dr.  Thus  the  helical  fields,  on  each  side  of  the  singular  surface, 
oppose  each  other.  As  the  fluid  inside  the  singular  surface  convects 
toward  one  side,  regions  of  plasma  with  opposing  helical  fluxes  are 
forced  together.  If  the  helical  flux  surfaces  maintain  their  topological 
structure,  as  they  would  in  ideal  MHD,  then  Rosenbluth,  Dagazian  and 
Rutherford  (Phys.  Fluids  1 6 ,  189^  (1973))  worked  out  a  nonlinear  theory 
for  the  assymptotic  displacement  of  the  plasma. •  They  did  this  by  balancing 
the  driving  force  against  the  additional  restraining  force  of  the  compressed 
field.  They  found  that  the  plasma  interior  is  displaced  only  slightly, 
for  instance  if  the  initial  constant  i|i  surfaces  are  shown  in  Fig.  (XIV  2a), 
the  final  constant  surfaces  might  appear  as  in  Fig.  (XIV  2b). 

Notice  however  that  in  this  final  state,  regions  of  opposite  helical 
field  (along  a-b)  are  forced  together.  If  reconnection  is  allowed  in 
this  region,  the  large  flux  compression  will  not  occur  and  the  plasmas 
sideways  displacement  can  proceed.  Let  us  look  into  this  reconnection 
more  carefully.  Since  the  magnetic  field  always  has  zero  divergence,  two 
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different  regions  can  reconnect  only  if  they  have  equal  helical  flux.  Thus, 
as  illustrated  in  Fig.  (XIV  la),  the  field  at  r^  can  connect  with  the  field 
at  r2«  If  the  layers  at  r^  and  r2  have  widths  dr.^  and  dr^,  the  fact  that 
the  fluxes  which  reconnect  are  equal  means 


Si*  -  n* 

fi  '» 


(XIV  3) 


These  two  flux  elements  will,  at  t  *  ",  come  to  rest  at  a  position  r  having 


width  dr.  Since  the  fluid  motion  is  incompressible 


rdr=  r,Jlr,  +  rt 


(XIV  k) 


Also,  since  the  reconnection  occurs  only  at  one  point,  (around  a-b),  the 


flux  is  undisturbed  in  most  of  the  flux  tube  so 


(XIV  5) 


Thus  as  the  field  reconnects,  the  total  helical  flux  in  the  reconnecting 
fluid  elements  is  conserved. 

It  is  not  difficult  to  follow  this  reconnection  process.  Initially 
flux  tubes  1  and  2  connect  as  shown  in  Fig.  (XIV  3  a  &  b)  and  an  island 
labeled  I  is  formed  opposite  to  the  direction  of  the  kink.  In  this  island, 
the  helical  field  is  clockwise,  as  it  is  outside  of  the  original  q  =  1 
surface.  Therefore  the  island  has  q  >  1.  Later  flux  elements  3  and  U 
reconnect  and  island  II  is  formed  outside  of  I  as  shown  in  Fig.  (XIV  3c). 
The  reconnection  process  completes  itself  when  flux  elements  5  and  6 
connect  to  form  island  III  (Fig.  (XIV  3d)  at  the  position  where 
i|i(r)  ■  <Ko)  Since  all  field  lines  are  now  clockwise,  q  is  everywhere 
greater  than  unity  and  the  plasma  has  returned  to  a  stable  state.  Thus 
the  flux  near  the  original  q=i  surface  has  moved  to  the  center  and  the 
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flux  initially  at  the  center  has  moved  to  the  outside.  The  <|/(r)  curve 
then  evolves  to  a  monotonically  decreasing  function  of  r  in  which  the 
area  between  two  values  of  ij;  is  the  same  at  t  =  »  as  at  t  =  0.  This  is 
illustrated  in  Fig.  (XIV  lb)  where  the  helical  flux  at  t  =  00  is  shown. 

Thus,  without  invoking  any  details  of  a  particular  m  =  n  =  1  instability, 
Kadomtsev  is  able  to  derive  the  field  configuration  of  a  final  stable 
state  from  an  initial  unstable  state.  Once  the  plasma  reaches  a  stable 
state,  the  channeling  of  the  current  into  the  hotter  regions  will  drive 
the  plasma  unstable,  as  discussed  in  Chapter  X.  Therefore  a  relaxation 
oscillation  will  ensue. 

It  is  worthwhile  to  briefly  compare  Kadomtsev's  theory  of  m  =  n  =  1 
instability  with  the  quasi-linear  theory  of  Chapter  X.  The  later  derived 
a  diffusion  equation  and  ultimately  flattened  the  current  profile  within 
the  q  =  1  surface.  Then  some  other  mechanism  was  invoked  to  spread  the 
current  beyond  the  original  q  =  1  surface  to  produce  a  plasma  with  q  >  1 
everywhere.  Kadomtsev's  theory  however  postulates  a  coupling  between  the 
plasma  inside  and  outside  of  the  q  -  1  surface  so  that  only  his  mechanism 
la  needed  to  produce  a  plasma  with  q  >  1  everywhere. 

We  now  turn  our  attention  to  a  very  brief  survey  of  numerical 
simulations  of  MHD  unstable  plasmas.  There  has  been  a  great  world  wide 
effort  in  such  simulations  recently  and  many  people  now  feel  that  this 
offers  the  best  hope  for  learning  about  the  nonlinear  theory  of  MHD 
instabilities.  Generally  speaking,  these  simulations  are  of  two  types. 

The  first  type,  pioneered  by  Rosenbluth  and  his  co-workers  at  Princeton, 
assumes  helical  symmetry  .incompressible  flow,  tokamak  ordering  and 
constant  density.  The  full  set  of  three  dimensional  MHD  equations  then 
reduce  to  two  equations  for  ifi  and  the  z  component  of  the  velocity  vector 


potential  in  the  two  dimensional  r,x  space.  The  second  approach  is  to  simply 
solve  the  full  set  of  MHD  equations  in  three  dimensions.  Here  the  plasma  is 
generally  assumed  to  fill  either  a  rectangular  solid,  or  else  a  torus  with 
rectangular  cross  section.  The  first  method  allows  much  greater  resolution, 
principally  because  the  helical  symmetry  reduces  the  dimensionality  from 
three  to  two.  The  second  method,  of  course,  has  much  greater  flexibility. 

The  m  =  1  internal  kink-tearing  mode  has  been  simulated  both  ways.  In 
each  case,  care  must  be  taken  to  allow  sufficient  resolution  in  the  singular 
layer.  This  has  necessitated  the  use  of  resistivity  much  larger  than  that 
existing  in  tokamak  plasmas. 

A  two  dimensional  simulation  of  the  kink  tearing  mode  was  performed  by 
Waddell,  Rosenbluth,  Monticello  and  White  (Nuclear  Fusion  l£,  528  (1978)). 

A  current  density 


Z 


[i*  NY 
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is  set  up  at  t  =  0.  J  is  chosen  so  that  q(r  =  0)  =  0.9  and  r  -  0.6a 

o  o 

where  a  is  the  radius  of  the  conducting  wall.  For  this  current  profile, 
the  radius  of  the  q  =  1  surface  is  at  r  =  0.2a  and  the  radius  at  which 
the  helical  flux  i|>  has  the  same  value  as  it  does  at  r  -  0  is  r^  «  0.3a. 
Island  formation  is  observed,  but  the  helical  flux  surfaces  become  quite 
complicated  very  quickly,  with  multiple  island  structure  developing.  What 
is  simpler  are  the  flow  patterns,  shown  at  four  times  in  Fig.  (XIV  h) . 

This  shows  that  the  basic  flow  pattern  of  an  m  =  1  mode  persists  well  into 
the  nonlinear  regime.  Kadomtsev's  theory  indicates  there  should  be 
reconnection  out  to  r^  =  0.3a.  The  perturbed  flow  and  perturbed  field 
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seem  to  be  limited  to  radii  somewhat  smaller  than  this,  but  accurate 


resolution  is  difficult.  Also,  Waddell  et  al  show  the  total  current  as  a 
function  of  radius,  and  this  is  shown  in  Fig.  (XIV  5).  Clearly,  the 
simple  quasi-linear  theory,  which  predicts  current  diffusion  within  the 
singular  surface,  is  reasonably  accurate. 

A  full  three  dimensional  simulation  of  this  instability  was  done  by 
Sykes  and  Wesson  (Phys.  Rev.  Lett,  37.,  1^0,  (1976)).  By  assuming  a 
resistivity  with  the  functional  form  n  ~  T  ”3/2,  they  were  able  to  simulate 
several  cycles  of  the  relaxation  oscillation.  In  Fig.  XIV  are  shown  various 
curves  of  constant  surfaces.  Here  it  can  be  seen  that  an  island  with 
q  >  1  grows  and  ultimately  displaces  the  original  flux  surfaces  having 
q  <  1.  However  different  cycles  of  the  relaxation  oscillation  are  not  all 
alike.  On  other  cycles,  the  original  island  may  return  and  displace  the 
newly  formed  island.  In  this  case  the  original  island  relaxes  to  a  state 
having  q  >  1. 

It  is  clear  that  num  erical  simulation  can  be  a  very  powerful  technique 
for  the  study  of  the  nonlinear  behavior  of  MHD  instabilities.  However,  as 
these  two  simulations  illustrate,  different  types  of  simulations  of  the 
same  process  do  not  necessarially  give  exactly  the  same  result.  Undoubtedly 
there  will  be  many  more  such  simulations  in  the  future.  In  fact,  many 
different  instabilities  have  already  been  studied  by  simulations,  including 
tearing  modes,  internal  kinks  and  free  boundary  kinks. 
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XIV  4)  The  time  evolution  of  the  flow  pattern  for  this  instability 


as  computed  by  Waddell  et  al 
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XIV  5)  The  time  evolution  of  the  toroidal  current  profile  as 
computed  by  Waddell  et  al. 


XIV  6)  The  time  evolution  of  the  magnetic  surfaces  as  computed 

by  Sykes  and  Wesson 


